library(dplyr)
library(tibble)
library(DBI)
library(sf)
library(gimms)
library(ncdf4)
library(rgdal)
library(stringr)
library(ggplot2)
library(RColorBrewer)
library(colorspace)
library(cowplot)
library(png)
library(spatialreg)
library(kableExtra)
library(reshape2)
library(grid)
library(gridExtra)
library(ggpubr)
source("f_plotspatial.R")
memory.limit(size = 160000)
## [1] Inf
library(captioner)
fig_nums <- captioner()
table_nums <- captioner(prefix = "Table")
The coastal drainage basins polygons were computed with a seamless digital elevation model for Fennoscandia (Finstad 2017). The model resulted in more than 200 000 coastal drainage basins, most of them being extremely small polygons covering the coastal regions. In order to simplify the model and reduce the amount of data, we keep only the 0,001 % of the coastal drainage basins that are the biggest. They account to more than 93% of the total surface covered by the initial amount of polygons.
con <- dbConnect(RPostgreSQL::PostgreSQL(),user = "camille.crapart", password = "camille",host = "vm-srv-wallace.vm.ntnu.no", dbname = "nofa")
query <- paste("SELECT gid,geom FROM ","\"Hydrography\"",".waterregions_dem_10m_nosefi",sep = "")
query.2 <- paste("SELECT * FROM ","\"Hydrography\"",".Fenoscandia_LakeRegister",sep = "")
water.regions <- st_read(con,query = query)
lake.register <- st_read(con,query = query.2)
saveRDS(water.regions,"water.regions.Rdata")
water.regions <- readRDS("water.regions.Rdata")
ggplot()+geom_sf(data=water.regions,aes(fill = gid))
water.regions <- readRDS("water.regions.Rdata")
water.regions$area <- st_area(water.regions) %>% as.numeric()
area1 <- sum(water.regions$area)
wr.quantile <- quantile(water.regions$area, probs = c(0.99,0.995,0.999))
wr.sf <- water.regions %>% filter(area > wr.quantile[3])
area2 <- sum(wr.sf$area)
ratio <- area2/area1*100
map.wr <- ggplot(wr.sf) + geom_sf(aes(fill=gid,col=gid))
ggsave(plot = map.wr, filename = "WR/wr.sf.png",dpi = "retina", width = 10, height = 15, units = "cm")
saveRDS(wr.sf,"WR/wr.sf.rds")
st_write(wr.sf, "wr.sf.shp")
NDVI values are extracted from the GIMMS NDVI3g dataset (The National Center for Atmospheric Research 2018), stored on http://poles.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, in the “ecocast” dataset and accessed via the “gimms” package (Detsch 2021) on Rstudio.
Data is taken bi-monthly (Tucker et al. 2005). Raster layers of all slices are downloaded using the gimmsdownload function, then monthly composite are calculated using the monthlyComposites function. The maximum NDVI value is taken for each pixel. Afterwards the NDVI values for each catchment polygon are extracted using the extract function from the raster package (Hijmans et al. 2018).
The values stored in the ecocast dataset are composed of the NDVI value and a flag value indicating the goodness of the data (Detsch 2021). All the NDVI values extracted, corresponding to the values for the studied catchments, had a flag of 1, indicating a good value. The NDVI value was retrieved from the stored value using the formula \(NDVI = floor(ndvi3g/10)/1000\) (Detsch 2021), and then the mean of the 3 summer values (June, July and August) was computed for each catchment.
ndvi.max <- readRDS("ndvi.max.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
summer.mean <- raster::stack(ndvi.max[[6]],ndvi.max[[7]],ndvi.max[[8]]) %>% mean()
summer.ndvi <- reclassify(summer.mean,cbind(-Inf, 0, NA), right=FALSE)
summer.scandinavia <- raster::crop(summer.ndvi,c(0,35,55,73))
writeRaster(summer.mean,"ndvi1994.tif",overwrite = T)
ndvi.wr <- raster::extract(summer.ndvi,wr.sf,fun = mean, df = T, sp = T, na.rm = T)
ndvi.wr$area <- NA
ndvi.wr$ndvi <- (floor(ndvi.wr$layer)/10)/1000
saveRDS(ndvi.wr, "WR/ndvi.wr.rds")
knitr::include_graphics("WR/wr.map.NDVI.png")
The runoff data was downloaded from the CORDEX database: https://portal.enes.org/data/enes-model-data/cordex. The data is available at https://esg-dn1.nsc.liu.se/projects/esgf-liu/.
The names of the variables in the CORDEX project follow the CF Metadata convention: http://cfconventions.org/
THe mean of the 30 years between 1970 and 2000 is used (p11 on EURO-CORDEX guidelines https://www.euro-cordex.net/imperia/md/content/csc/cordex/guidance_for_euro-cordex_climate_projections_data_use__2021-02_1_.pdf). This time period matches the temperature and precipitation data provided by WorldClim.
runoff.raster <- "CORDEX/mrros_Lmon_CNRM-CM6-1-HR_historical_r1i1p1f2_gr_185001-201412.nc"
runoff.nc <- nc_open(runoff.raster)
runoff.stack <- raster::stack(runoff.raster, bands = c(1441:1812))
runoff.mean <- runoff.stack %>% mean()
wr.sf <- readRDS("WR/wr.sf.rds")
runoff.wr95 <- raster::extract(runoff.mean, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(runoff.wr95) <- c("gid","area","runoff")
saveRDS(runoff.wr95,"WR/runoff.wr95.rds")
knitr::include_graphics("WR/wr.map.Runoff.png")
The Land Cover data was downloaded from https://land.copernicus.eu/pan-european/corine-land-cover/lcc-2000-2006. The previous version (1995-2000) excludes Norway so the 2000 version was preferred.
The land cover data was downloaded as a raster (.tiff file), which included the legend recording the code for each land use. The raster and legend were assembled in QGIS before being imported as a shapefile in R. The codes (corresponding to the line in the extracted dataframe) of the categories of interest were:
Arable land:
Bogs:
Forest:
Bare:
wr.sf <- readRDS("WR/wr.sf.rds")
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
wr.corine <- st_transform(wr.sf, crs(corine))
extract_corine <- function(x){
clc <- extract(corine,x)
saveRDS(clc,file = paste("WR/corine/clc",x$gid,"rds",sep="."))
}
wr.list <- split(wr.corine,seq(nrow(wr.corine)))
lapply(wr.list[1143:1440],extract_corine)
#clc.test <- readRDS("WR/clc/clc.16711649.rds")
clc.files <- list.files(path = "WR/corine/", pattern = "clc.*.rds", full.names = T)
list.clc <- sapply(clc.files,readRDS)
list.clc.wr <- sapply(list.clc, as.integer)
names.list <- str_extract(clc.files,"\\d+")
names(list.clc.wr) <- names.list
saveRDS(list.clc.wr,"WR/corine/list.clc.wr.rds")
list.clc.wr <- readRDS("WR/corine/list.clc.wr.rds")
clc.tab.wr <- sapply(list.clc.wr, function(x) tabulate(x,45)) %>% t()
clc.tab.area.wr <- clc.tab.wr*prod(res(corine))
catchment.area <- rowSums(clc.tab.area.wr)
clc.tab.prop.wr <- sweep(clc.tab.area.wr,1,catchment.area, FUN = "/")
saveRDS(clc.tab.prop.wr,"WR/corine/clc.tab.prop.wr.rds")
clc.tab.prop.wr <- readRDS("WR/corine/clc.tab.prop.wr.rds") %>% as.data.frame() %>% tibble::rownames_to_column(var = "gid")
bogs <- clc.tab.prop.wr %>% dplyr::select(c("gid","V36")) %>% setNames(c("gid","bogs"))
saveRDS(bogs,"WR/corine/bogs.wr.rds")
arable <- data.frame( gid = clc.tab.prop.wr$gid, arable = rowSums(select(clc.tab.prop.wr, c("V12","V16","V18","V19","V20"))))
saveRDS(arable,"WR/corine/arable.wr.rds")
forest <- data.frame( gid = clc.tab.prop.wr$gid, forest = rowSums(select(clc.tab.prop.wr, c("V23","V24","V25"))))
saveRDS(forest,"WR/corine/forest.wr.rds")
bare <- data.frame( gid = clc.tab.prop.wr$gid, bare = rowSums(select(clc.tab.prop.wr, c("V31","V32","V33"))))
saveRDS(bare,"WR/corine/bare.wr.rds")
water <- data.frame( gid = clc.tab.prop.wr$gid, water = rowSums(select(clc.tab.prop.wr, c("V40","V41","V42","V43","V44"))))
saveRDS(water,"WR/corine/water.wr.rds")
list.grob <- c("WR/wr.map.Bog.png", "WR/wr.map.Arable.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 1, ncol = 2, labels = c("a)","b)"), label_size = 10)
Temperature and precipitation are downloaded from the WorldClim database: https://worldclim.org/data/index.html
Temperature in historical Worldclim is in °C * 10 (to reduce file size). Precipitation are average monthly precipitation in mm. Here we use the mean monthly precipitation and temperature for the period 1970-2000.
wr.sf <- readRDS("WR/wr.sf.rds")
bio.fennoscandia <- raster::getData(name = "worldclim", var = "bio", res = 2.5)
temp.wr <- raster::extract(bio.fennoscandia$bio1, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(temp.wr) <- c("gid","area","temp")
temp.wr$temp <- temp.wr$temp / 10
saveRDS(temp.wr,"WR/temp.wr.rds")
#temp.wr <- readRDS("WR/temp.wr.rds")
prec.wr <- raster::extract(bio.fennoscandia$bio12, wr.sf, fun = mean, df = T, sp = T, na.rm = T)
names(prec.wr) <- c("gid","area","prec")
saveRDS(prec.wr,"WR/prec.wr.rds")
#prec.wr <- readRDS("WR/prec.wr.rds")
tp <- cbind(wr.sf[,1],temp.wr[,3]) %>% cbind(prec.wr[,3])
saveRDS(tp,"WR/tp.rds")
list.grob <- c("WR/wr.map.Temp.png", "WR/wr.map.Precip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow= 1, ncol = 2, labels = c("a)", "b)"), label_size = 10)
The nitrogen deposition data was extracted from the EMEP database (https://emep.int/mscw/mscw_moddata.html#Comp). The data from 2000 was used as it is the first model available.
N-deposition (NOx and NH3) is the sum of:
library(ncdf4)
EMEP_file <- "EMEP/EMEP01_rv4.42_year.2000met_2000emis_rep2021.nc"
wr.sf <- readRDS("WR/wr.sf.rds")
woxn <- raster(EMEP_file, varname = "WDEP_OXN")
doxn <- raster(EMEP_file, varname = "DDEP_OXN_m2Grid")
wrdn <- raster(EMEP_file, varname = "WDEP_RDN")
drdn <- raster(EMEP_file, varname = "DDEP_RDN_m2Grid")
woxn_wr <- extract(woxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(woxn_wr) <- c("gid","area","woxn")
saveRDS(woxn_wr,"woxn_wr.rds")
doxn_wr <- extract(doxn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(doxn_wr) <- c("gid","area","doxn")
saveRDS(doxn_wr,"doxn_wr.rds")
wrdn_wr <- extract(wrdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(wrdn_wr) <- c("gid","area","wrdn")
saveRDS(wrdn_wr,"wrdn_wr.rds")
drdn_wr <- extract(drdn,wr.sf, sp = T, df = T, na.rm = T, fun = mean)
names(drdn_wr) <- c("gid","area","drdn")
saveRDS(drdn_wr,"drdn_wr.rds")
ndep.wr <- merge(woxn_wr,doxn_wr, by = "gid") %>% merge(wrdn_wr, by = "gid") %>% merge(drdn_wr, by = "gid")
ndep.wr$tndep <- rowSums(ndep.wr@data[,c(3:6)])
saveRDS(ndep.wr,"WR/ndep.wr.rds")
knitr::include_graphics("WR/wr.map.TNdep.png")
The altitude data was accessed through the “getData” function, and comes from the NASA digital elevation model (SRTM 90 m). The slope was then computed using the “terrain” function from the raster package.
wr.sf <- readRDS("WR/wr.sf.rds")
alt.nor <- raster::getData("alt",country = "NOR", mask = T)
alt.swe <- raster::getData("alt",country = "SWE", mask = T)
alt.fin <- raster::getData("alt",country = "FIN", mask = T)
alt.list <- list(alt.nor,alt.swe,alt.fin)
alt <- do.call(raster::merge,alt.list)
alt.wr <- raster::extract(alt, dplyr::select(wr.sf, c("gid","geom")), sp = T, fun = mean, df = T, na.rm = T)
names(alt.wr) <- c("gid","alt")
saveRDS(alt.wr,"WR/alt.wr.rds")
slope.nor <- terrain(alt.nor, unit = "degrees")
slope.swe <- terrain(alt.swe, unit = "degrees")
slope.fin <- terrain(alt.fin, unit = "degrees")
slope.list <- list(slope.nor,slope.swe,slope.fin)
slope <- do.call(raster::merge,slope.list)
slope.nona <- reclassify(slope,cbind(NA,100), right = NA)
slope.wr <- raster::extract(slope,wr.sf, sp = T, df = T, fun = mean, na.rm = T)
names(slope.wr) <- c("gid","area","slope")
saveRDS(slope.wr, "WR/slope.wr.rds")
The data was gathered and merged into a single dataframe, based on the catchment identifyer (“gid”).
ndvi.wr <- readRDS("WR/ndvi.wr.rds")
ndvi.wr$area <- NULL
runoff.wr95 <- readRDS("WR/runoff.wr95.rds")
runoff.wr95$area <- NULL
bogs.wr <- readRDS("WR/corine/bogs.wr.rds")
forest.wr <- readRDS("WR/corine/forest.wr.rds")
arable.wr <- readRDS("WR/corine/arable.wr.rds")
bare.wr <- readRDS("WR/corine/bare.wr.rds")
water.wr <- readRDS("WR/corine/water.wr.rds")
#clc.wr.sum <- readRDS("WR/clc/clc.wr.sum.Rdata")
alt.wr <- readRDS("WR/alt.wr.rds")
slope.wr <- readRDS("WR/slope.wr.rds")
slope.wr$area <- NULL
tp <- readRDS("WR/tp.rds")
ndep.wr <- readRDS("WR/ndep.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
all.wr.sf <- merge(wr.sf,runoff.wr95,by = "gid") %>% merge(ndvi.wr,by="gid") %>% merge(bogs.wr,by = "gid") %>% merge(forest.wr, by = "gid") %>% merge(arable.wr, by = "gid") %>% merge(bare.wr, by = "gid") %>% merge(water.wr, by = "gid") %>% merge(st_drop_geometry(tp),by = "gid") %>% merge(alt.wr, by = "gid") %>% merge(slope.wr, by = "gid") %>% merge(ndep.wr, by = "gid")
saveRDS(all.wr.sf, "WR/all.wr.sf.rds")
The data was then cleaned and the centroid of each coastal drainage basin was computed.
all.wr.sf <- readRDS("WR/all.wr.sf.rds")
wr.sf.95 <- all.wr.sf %>% dplyr::filter(ndvi != 0) %>% filter(runoff > 0) %>% filter(is.na(bogs) == F) %>% filter(is.na(alt) == F) %>% filter(is.na(slope) == F)
wr.sf.95 <- wr.sf.95[!duplicated(wr.sf.95$gid),]
wr.sf.95$Runoff <- wr.sf.95$runoff *365*24*60*60 # from kg/m2/s to kg/m2/year or mm/y
wr.sf.95$logRunoff <- log(wr.sf.95$Runoff)
sfc_point_to_matrix = function(x) {
matrix(
unlist(x, use.names = FALSE),
nrow = length(x),
byrow = TRUE,
dimnames = list(1:length(x))
)
}
wr.centroid <- st_centroid(wr.sf.95, of.largest.polygon = T)
coord <- sfc_point_to_matrix(wr.centroid$geometry) %>% as.data.frame() %>% setNames(c("Longitude","Latitude"))
wr.sf.95 <- cbind(wr.sf.95, coord)
wr.sf.95$area <- st_area(wr.sf.95)
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")
noforest <- filter(wr.sf.95,forest == 0) #%>% dplyr::select(noforest,c("gid","area"))
#ggplot(noforest[which.max(noforest$area),])+geom_sf(aes(fill=as.factor(gid)))
weird.gid <- noforest$gid[which.max(noforest$area)]
weird.cat <- wr.sf.95 %>% filter(gid == weird.gid) %>% st_transform(CRS("EPSG:3035"))
corine <- raster::raster("Corine_Land_Cover/U2006_CLC2000_V2020_20u1.tif")
#weird.clc <- raster::extract(corine,weird.cat)
#saveRDS(weird.clc,"WR/clc/weird.clc.rds")
weird.clc <- readRDS("WR/clc/weird.clc.rds")
weird.tab <- sapply(weird.clc, function(x) tabulate(x,45))
weird.tabl.area <- weird.tab*prod(res(corine))
catchment.area <- colSums(weird.tabl.area)
weird.tab.prop <- sweep(weird.tabl.area,2,catchment.area, FUN = "/")
wr.sf.95$bogs[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[36]
wr.sf.95$forest[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[24]+weird.tab.prop[24]+weird.tab.prop[25]
wr.sf.95$arable[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[12]+weird.tab.prop[16]+weird.tab.prop[18]+weird.tab.prop[19]+weird.tab.prop[20]
wr.sf.95$glacier[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[34]
wr.sf.95$bare[which(wr.sf.95$gid == weird.cat$gid)] <- weird.tab.prop[31]+weird.tab.prop[32]+weird.tab.prop[33]
saveRDS(wr.sf.95,"WR/wr.sf.95.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
toreplace <- which(names(wr.sf.95) %in% c("ndvi","bogs","forest","arable","bare","water","temp","prec","alt","slope","tndep"))
names(wr.sf.95)[toreplace] <- c("NDVI","Bog","Forest","Arable","Bare","Water","Temp","Precip","Alt","Slope","TNdep")
wr.sf.95$area.x <- NULL
wr.sf.95$area.y <- NULL
saveRDS(wr.sf.95, "WR/wr.sf.95.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
m <- c("NDVI","Runoff","Precip","Temp","Bog","Forest","Arable","TNdep")
u <- c("NDVI on scale 0 to 1", "Mean Runoff, mm/y","Mean monthly precipitation,\n mm","Mean annual temperature,\n °C","Proportion of bogs","Proportion of Forest","Proportion of\narable land","Total Nitrogen Deposition (mg/m2)")
couleurs <- c("Greens","Blues","Blues","RdYlBu","YlOrBr","RdYlGn","RdYlGn","PuRd")
dir <- c(T,F,F,T,T,F,T,F)
for(i in m){
mysf <- wr.sf.95
map <- ggplot(mysf)+geom_sf(aes(fill = .data[[i]], col = .data[[i]]))+
scale_fill_continuous_divergingx(palette = "RdYlBu",
aesthetics = c("colour","fill"),
rev = dir[grep(i,m)],
name = paste(u[grep(i,m)],"\n(",i,")",sep=""))+
theme_minimal(base_size = 9)+
theme(legend.position = "bottom")
filename_i <- paste("WR/wr.map",i,"png",sep = ".")
ggsave(plot = map, filename = filename_i, dpi = "retina", width = 10, height = 12, units = "cm")
}
wr.sf.95 <- readRDS("WR/wr.sf.95.rds") %>% st_drop_geometry()
wr.sf.95$area_km2 <- wr.sf.95$area * 10^(-6)
wr.min <- lapply(wr.sf.95, min) %>% as.data.frame()
wr.max <- lapply(wr.sf.95, max) %>% as.data.frame()
wr.median <- lapply(wr.sf.95, median) %>% as.data.frame()
wr.mean <- lapply(wr.sf.95, mean) %>% as.data.frame()
wr.sd <- lapply(wr.sf.95, sd) %>% as.data.frame()
wr.sum <- rbind(wr.min, wr.max, wr.median, wr.mean, wr.sd)
rownames(wr.sum) <- c("Min","Max","Median","Mean","Sd")
wr.sum %>% dplyr::select(c("NDVI", "Forest", "Arable","Temp","Precip","TNdep","Runoff","area_km2")) %>% kable(digits = 2) %>%
kable_styling(bootstrap_options = "bordered")
| NDVI | Forest | Arable | Temp | Precip | TNdep | Runoff | area_km2 | |
|---|---|---|---|---|---|---|---|---|
| Min | 0.11 | 0.00 | 0.00 | -4.24 | 434.55 | 63.88 | 55.17 | 27.3728 [m^2] |
| Max | 0.90 | 0.96 | 0.92 | 8.13 | 2842.75 | 2387.52 | 971.21 | 49832.2278 [m^2] |
| Median | 0.73 | 0.38 | 0.01 | 3.24 | 671.00 | 432.82 | 250.68 | 76.2258 [m^2] |
| Mean | 0.68 | 0.39 | 0.08 | 3.05 | 966.33 | 586.44 | 339.81 | 751.3046 [m^2] |
| Sd | 0.15 | 0.25 | 0.15 | 2.97 | 582.82 | 491.72 | 222.63 | 3752.6959 [m^2] |
We model the mean TOC concentration in all the water region in 1995 based on the simple SEM model of the previous section, using NDVI, bogs and runoff as predictors.
The neighbor matrix of the coastal drainage basins polygons were computed using their centroids.
wr.sf.95 <- readRDS("wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T)
wr.neigh.set <- knearneigh(wr.centroid, k = 100)
wr.neigh.nb <- knn2nb(wr.neigh.set)
wr.kmat <- nb2listw(wr.neigh.nb)
saveRDS(wr.kmat,"WR/wr.kmat.rds")
The TOC concentration in coastal drainage basins was modeled using the Spatial Error Linear Model with 5 predictors:
library(spatialreg)
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred <- predict(sem.fennoscandia.5,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
sem.pred.wr <- cbind(wr.sf.95,wr.pred$fit)
sem.pred.wr$TOC95 <- exp(sem.pred.wr$wr.pred.fit)
map.fit <- ggplot()+geom_sf(data = st_as_sf(sem.pred.wr), aes(fill = TOC95, col = TOC95))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled TOC 1995 (mg/L)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.fit, filename = "WR/wr.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
map.log <- ggplot()+geom_sf(data = sem.pred.wr, aes(fill = wr.pred.fit, col = wr.pred.fit))+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1,
name = "Modelled log(TOC)")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
ggsave(plot = map.log, filename = "WR/wr.log.toc.png",dpi = "retina", width = 6, height = 8, units = "cm")
saveRDS(sem.pred.wr,"WR/sem.pred.wr.rds")
knitr::include_graphics("WR/wr.log.toc.png")
The forecast of the temperature and precipitation in the periods 2041-2060 and 2081-2100, for SSP 1-2.6 and SSP 3-7.0, were extracted based on the Future Worldclim dataset.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 1)
ssp370.cnrm.temp <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 1)
ssp126.cnrm.temp.2050.crop <- raster::crop(ssp126.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2050.png")
plot(ssp126.cnrm.temp.2050.crop)
dev.off()
ssp370.cnrm.temp.2050.crop <- raster::crop(ssp370.cnrm.temp,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2050.png")
plot(ssp370.cnrm.temp.2050.crop)
dev.off()
ssp126.temp.wr <- raster::extract(ssp126.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.wr) <- c("gid", "temp126")
saveRDS(ssp126.temp.wr,"WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- raster::extract(ssp370.cnrm.temp,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.wr) <- c("gid", "temp370")
saveRDS(ssp370.temp.wr,"WR/ssp370.temp.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 1)
raster.ssp370.temp.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 1)
ssp126.cnrm.temp.2100.crop <- raster::crop(raster.ssp126.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp126.temp.2100.png")
plot(ssp126.cnrm.temp.2100.crop)
dev.off()
ssp370.cnrm.temp.2100.crop <- raster::crop(raster.ssp370.temp.2100,c(0,35,55,73))
png("worldclim_future/ssp370.temp.2100.png")
plot(ssp370.cnrm.temp.2100.crop)
dev.off()
ssp126.temp.2100 <- raster::extract(raster.ssp126.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm =T)
names(ssp126.temp.2100) <- c("gid", "temp126.2100")
saveRDS(ssp126.temp.2100,"WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- raster::extract(raster.ssp370.temp.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.temp.2100) <- c("gid", "temp370.2100")
saveRDS(ssp370.temp.2100,"WR/ssp370.temp.2100.rds")
listgrob <- list("worldclim_future/ssp126.temp.2050.png", "worldclim_future/ssp370.temp.2050.png","worldclim_future/ssp126.temp.2100.png","worldclim_future/ssp370.temp.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp126_2041-2060.tif", band = 12)
ssp370.cnrm.prec <- raster::raster("worldclim_future/wc2.1_2.5m_bioc_CNRM-CM6-1_ssp370_2041-2060.tif", band = 12)
ssp126.cnrm.prec.crop <- raster::crop(ssp126.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp126.2050.prec.png")
plot(ssp126.cnrm.prec.crop)
dev.off()
ssp370.cnrm.prec.crop <- raster::crop(ssp370.cnrm.prec,c(0,35,55,73))
png("worldclim_future/ssp370.2050.prec.png")
plot(ssp370.cnrm.prec.crop)
dev.off()
ssp126.prec.wr <- raster::extract(ssp126.cnrm.prec,select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.wr) <- c("gid", "prec126")
saveRDS(ssp126.prec.wr,"WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- raster::extract(ssp370.cnrm.prec,dplyr::select(wr.sf.95,c("gid","geometry")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.wr) <- c("gid", "prec370")
saveRDS(ssp370.prec.wr,"WR/ssp370.prec.wr.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp126_2081-2100.tif", band = 12)
raster.ssp370.prec.2100 <- raster::raster("worldclim_future/wc2.1_30s_bioc_CNRM-CM6-1-HR_ssp370_2081-2100.tif", band = 12)
ssp126.cnrm.prec.2100.crop <- raster::crop(raster.ssp126.prec.2100,c(0,35,55,73))
png("worldclim_future/ssp126.2100.prec.png")
plot(ssp126.cnrm.prec.2100.crop)
dev.off()
png("worldclim_future/ssp370.2100.prec.png")
ssp370.cnrm.prec.2100.crop <- raster::crop(raster.ssp370.prec.2100,c(0,35,55,73))
plot(ssp370.cnrm.prec.2100.crop)
ssp126.prec.2100 <- raster::extract(raster.ssp126.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp126.prec.2100) <- c("gid", "prec126.2100")
saveRDS(ssp126.prec.2100,"WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- raster::extract(raster.ssp370.prec.2100,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T, na.rm = T)
names(ssp370.prec.2100) <- c("gid", "prec370.2100")
saveRDS(ssp370.prec.2100,"WR/ssp370.prec.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
ssp126.temp.wr <- readRDS("WR/ssp126.temp.wr.rds")
ssp370.temp.wr <- readRDS("WR/ssp370.temp.wr.rds")
ssp126.prec.wr <- readRDS("WR/ssp126.prec.wr.rds")
ssp370.prec.wr <- readRDS("WR/ssp370.prec.wr.rds")
ssp126.temp.2100 <- readRDS("WR/ssp126.temp.2100.rds")
ssp370.temp.2100 <- readRDS("WR/ssp370.temp.2100.rds")
ssp126.prec.2100 <- readRDS("WR/ssp126.prec.2100.rds")
ssp370.prec.2100 <- readRDS("WR/ssp370.prec.2100.rds")
# SSP126
df.ssp126 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.wr, by = "gid") %>% merge(ssp126.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126$logPrecip <- df.ssp126$Precip %>% log()
saveRDS(df.ssp126, "WR/df.ssp126.rds")
# SSP 370
df.ssp370 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.wr, by = "gid") %>% merge(ssp370.prec.wr, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370$logPrecip <- df.ssp370$Precip %>% log()
saveRDS(df.ssp370, "WR/df.ssp370.rds")
# SSP126 2100
df.ssp126.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp126.temp.2100, by = "gid") %>% merge(ssp126.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp126.2100$logPrecip <- df.ssp126.2100$Precip %>% log()
saveRDS(df.ssp126.2100, "WR/df.ssp126.2100.rds")
# SSP 370 2100
df.ssp370.2100 <- dplyr::select(wr.sf.95, c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff")) %>% merge(ssp370.temp.2100, by = "gid") %>% merge(ssp370.prec.2100, by = "gid") %>% st_drop_geometry() %>% setNames(c("gid","area","Alt","Slope","Latitude","Bare", "Arable", "Runoff","Temp","Precip"))
df.ssp370.2100$logPrecip <- df.ssp370.2100$Precip %>% log()
saveRDS(df.ssp370.2100, "WR/df.ssp370.2100.rds")
listgrob <- list("worldclim_future/ssp126.2050.prec.png", "worldclim_future/ssp370.2050.prec.png","worldclim_future/ssp126.2100.prec.png","worldclim_future/ssp370.2100.prec.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("a) SSP126-2050","b) SSP370-2050","c) SSP126-2100","d) SSP370-2100"), label_size = 10)
We model NDVI based on temperature and precipitation, and use this model to predict the future NDVI in 2050. Choice of future data: ssp 126 as best scenario, ssp 370 as worse scenario (Hausfather 2019; CORDEX 2021). CRNM according to Raju and Kumar (2020).
Using a SELM model for NDVI, we notice that high residuals are observed in catchments where original NDVI is lower than 0,4. Myers-Smith et al. (2020) show that the relationship between NDVI and the Arctic biomass is not linear: under 0,4, an increase in NDVI reflects the vegetalization of bare ground, while over 0,4, an increase of NDVI reflects an increased productivity in established forests.
Several models were tested to model the NDVI: a linear model, a spatial error linear model, a polynomial model, and a polynomial model combined with a binomial model. The NDVI was fitted on the 1995 data from the Northern Lakes Survey and tested on on the coastal drainage basins.
library(betareg)
library(mgcv)
library(VGAM)
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
fennoscandia$Latitude <- fennoscandia$latitude
fen.kmat <- readRDS("k.neigh.w.rds")
ndvi.cor <- cor(dplyr::select(st_drop_geometry(fennoscandia), c("NDVI", "Temp", "Precip", "Arable", "Bare", "Latitude", "Slope","Runoff")))
corrplot::corrplot(ndvi.cor, method = "number")
par(mfrow = c(2,3))
plot(NDVI~Temp+log(Precip)+Arable+Bare+Latitude+Slope, data = st_drop_geometry(fennoscandia))
par(mfrow = c(3,2))
#LM
summary(m1 <- lm(NDVI ~ Temp + log(Precip), data = fennoscandia))
summary(m1 <- step(m1))
plot(predict(m1), fennoscandia$NDVI, pch = ".", main="a ) NDVI - LM", ylab="Observed NDVI", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.4, label = paste("r2 =", round(summary(m1)$adj.r.squared,2), sep= ""))
#SELM
fen.kmat <- readRDS("k.neigh.w.rds")
selm.model <- errorsarlm(formula = NDVI~Temp+logPrecip, data = fennoscandia, listw = fen.kmat)
saveRDS(selm.model, "WR/selm.ndvi.rds")
selm.model <- readRDS("WR/selm.ndvi.rds")
AIC <- 2*(3-1)-2*ndvi.model$LL
plot(ndvi.model$fitted.value, fennoscandia$NDVI, pch = ".", main="b) NDVI - SELM", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.8, y = 0.5, label = paste("r2 =", round(cor(ndvi.model$fitted.values, fennoscandia$NDVI)^2,2), sep= ""))
# quadratic
summary(m2 <- lm(NDVI ~ (Temp + I(Temp^2)) * (log(Precip) + I(log(Precip)^2)), data = fennoscandia))
summary(m2 <- step(m2))
plot(predict(m2), fennoscandia$NDVI, pch = ".", main="c) NDVI - quadratic model", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m2)$adj.r.squared,2), sep= ""))
# poly
summary(m3 <- lm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia))
summary(m3 <- step(m3))
plot(predict(m3), fennoscandia$NDVI, pch = ".", main="d) NDVI - polynomial ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r2 =", round(summary(m3)$adj.r.squared,2), sep= ""))
saveRDS(m3,"WR/poly.ndvi.rds")
# glm
summary(m4 <- glm(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, family = "quasibinomial"))
#summary(m4 <- step(m4))
plot(predict(m4, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
text(x = 0.7, y = 0.4, label = paste("r =", cor(fennoscandia$NDVI, m4$fitted.values), sep= ""))
saveRDS(m4, "WR/glm.ndvi.rds")
#betareg
summary(m5 <- betareg(NDVI ~ poly(cbind(Temp, logPrecip), degree = 2), data = fennoscandia, link = "logit"))
plot(predict(m5, type = "response"), fennoscandia$NDVI, pch = ".", main="e) NDVI - glm ", ylab="Observed", xlab="Predicted",abline(0,1))
saveRDS(m5, "WR/betareg.ndvi.rds")
The SELM and the polynomial model had the best results. We tested them on the coastal drainage basins dataset.
k.neigh.w <- readRDS("k.neigh.w.rds")
fennoscandia <- readRDS("fennoscandia.rds")
fennoscandia.33 <- readRDS("fennoscandia.33.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
selm.ndvi.fen <- predict(selm.ndvi,wr.sf.95,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.selm.ndvi.model <- wr.sf.95 %>% st_as_sf() %>% mutate(ndvi.fit = selm.ndvi.fen$fit, ndvi.residuals = wr.sf.95$NDVI - selm.ndvi.fen$fit)
wr.selm.ndvi.model$area <- st_area(wr.selm.ndvi.model) %>% as.numeric()
ggplot(wr.selm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.selm.ndvi.model$NDVI,wr.selm.ndvi.model$ndvi.fit)^2,2)), x = 0.25, y = 0.5, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for SELM")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.selm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals.selm <- ggplot()+
geom_sf(data = wr.selm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "SELM residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for SELM")
ggsave(plot = wr.residuals.selm, filename = "WR/wr.selm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
fennoscandia <- readRDS("fennoscandia.33.rds")
poly.ndvi.wr <- predict(poly.ndvi,wr.sf.95, type = "response")
wr.poly.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = poly.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - poly.ndvi.wr)
ggplot(wr.poly.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.poly.ndvi.model$NDVI,wr.poly.ndvi.model$ndvi.fit),2)), x = 0.75, y = 0.4, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for polynomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.poly.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.poly.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "Polynomial residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Polynomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.poly.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
glm.ndvi.wr <- predict(glm.ndvi,wr.sf.95, type = "response")
wr.glm.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = glm.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - glm.ndvi.wr)
ggplot(wr.glm.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.glm.ndvi.model$NDVI,wr.glm.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.glm.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.glm.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "glm residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.glm.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- log(wr.sf.95$Precip)
betareg.ndvi.wr <- predict(betareg.ndvi,wr.sf.95, type = "response")
wr.betareg.ndvi.model <- wr.sf.95 %>% mutate(ndvi.fit = betareg.ndvi.wr, ndvi.residuals = wr.sf.95$NDVI - betareg.ndvi.wr)
ggplot(wr.betareg.ndvi.model) + geom_point(aes(y = ndvi.fit, x = NDVI, col = Latitude), size = 1)+
scale_color_distiller(type = "div", palette = "RdYlBu", direction = 1, aesthetics = c("col"))+
labs(x = "NDVI", y = "Predicted NDVI")+
annotate(geom = "label", label = paste("r = ", round(cor(wr.betareg.ndvi.model$NDVI,wr.betareg.ndvi.model$ndvi.fit),2)), x = 0.8, y = 0.3, size = 2)+
geom_abline(slope = 1, intercept = 0, col = "red")+
labs(title = "Observed vs predicted NDVI for binomial model")+
theme_minimal(base_size = 6)
ggsave(filename = "WR/wr.ndvi.betareg.results.png", dpi = "retina", width = 10, height = 8, units = "cm")
wr.residuals <- ggplot()+
geom_sf(data = wr.betareg.ndvi.model, aes(fill = ndvi.residuals, col = ndvi.residuals))+
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, aesthetics = c("fill", "col"), name = "betareg residuals")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")+
labs(title = "Map of residuals for Binomial model")
ggsave(plot = wr.residuals, filename = "WR/wr.betareg.ndvi.residuals.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob1 <- list("WR/wr.ndvi.selm.results.png","WR/wr.ndvi.poly.results.png","WR/wr.ndvi.glm.results.png","WR/wr.ndvi.betareg.results.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob1, ncol = 2, nrow = 2, labels = c("a)","b)","c)","d)"))
listgrob2 <- list("WR/wr.selm.ndvi.residuals.png","WR/wr.poly.ndvi.residuals.png", "WR/wr.glm.ndvi.residuals.png","WR/wr.betareg.ndvi.residuals.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob2, ncol= 2, nrow = 2,labels = c("e)","f)","g)","h)"))
selm.ndvi <- readRDS("WR/selm.ndvi.rds")
poly.ndvi <- readRDS("WR/poly.ndvi.rds")
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
betareg.ndvi <- readRDS("WR/betareg.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
xtemp <- seq(min(fennoscandia$Temp), max(fennoscandia$Temp), 0.1)
yprecip <- seq(min(fennoscandia$logPrecip), max(fennoscandia$logPrecip), 0.1)
data.fit <- expand.grid(Temp = xtemp, logPrecip = yprecip)
mtrx3d <- predict(selm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp1 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot SELM")+
theme_minimal(base_size = 10)
mtrx3d <- predict(poly.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp2 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot polynomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(glm.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp3 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot binomial model")+
theme_minimal(base_size = 10)
mtrx3d <- predict(betareg.ndvi, newdata = data.fit, type = "response")
mtrx.melt <- cbind(data.fit, mtrx3d)
names(mtrx.melt) <- c("Temp","logPrecip","NDVI")
cp4 <- ggplot(mtrx.melt,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..), col = "black")+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = fennoscandia, aes(x=Temp, y = logPrecip, col = NDVI), size = 1)+
scale_fill_viridis_c(aesthetics = c("fill", "col"), name = "NDVI")+
labs(title = "Contour plot beta model")+
theme_minimal(base_size = 10)
grid.arrange(cp1,cp2,cp3,cp4, nrow = 2, ncol = 2)
The binomial model preforms better: low residuals and prediction range between 0 and 1. This model will be used to forecast NDVI in 2050 and in 2100.
df.ssp126 <- readRDS("WR/df.ssp126.rds")
df.ssp370 <- readRDS("WR/df.ssp370.rds")
df.ssp126.2100 <- readRDS("WR/df.ssp126.2100.rds")
df.ssp370.2100 <- readRDS("WR/df.ssp370.2100.rds")
ndvi.model <- readRDS("WR/betareg.ndvi.rds")
wr.kmat <- readRDS("WR/wr.kmat.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.sf.95$logPrecip <- wr.sf.95$Precip %>% log()
#NDVI 126 2050
# ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.fit <- predict(object = ndvi.model, newdata = df.ssp126, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126 <- cbind(df.ssp126, ndvi.ssp126.fit)
saveRDS(ndvi.ssp126,"WR/ndvi.ssp126.rds")
# NDVI 370 2050
# ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.fit <- predict(object = ndvi.model, newdata = df.ssp370, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370 <- cbind(df.ssp370, ndvi.ssp370.fit)
saveRDS(ndvi.ssp370,"WR/ndvi.ssp370.rds")
#NDVI 126 2010
# ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% as.data.frame() %>% setNames(c("ndvi","trend","signal"))
#
# ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp126.2100.fit <- predict(object = ndvi.model, newdata = df.ssp126.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp126.2100 <- cbind(df.ssp126.2100, ndvi.ssp126.2100.fit)
saveRDS(ndvi.ssp126.2100,"WR/ndvi.ssp126.2100.rds")
#NDVI 370 2100
# ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, listw = wr.kmat, pred.type = "TS", legacy.mixed = F, zero.policy = T) %>% setNames(c("ndvi","trend","signal"))
# ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit) %>% setNames(c("gid","area","Alt","Slope","Latitude","Temp","Precip","NDVI","trend","signal","geometry"))
ndvi.ssp370.2100.fit <- predict(object = ndvi.model, newdata = df.ssp370.2100, type = "response") %>% as.data.frame() %>% setNames("NDVI")
ndvi.ssp370.2100 <- cbind(df.ssp370.2100, ndvi.ssp370.2100.fit)
saveRDS(ndvi.ssp370.2100,"WR/ndvi.ssp370.2100.rds")
In the historical dataset in Worldclim, the runoff is given as a yearly average in kg/m2/s, or mm/s.
In the CNRM-CM6 forecast we get the monthly forecast for runoff (in mm/s) and convert it in mm/y.
listgrob <- list("CORDEX_future/runoff.ssp126.2050.png", "CORDEX_future/runoff.ssp370.2050.png","CORDEX_future/runoff.ssp126.2100.png","CORDEX_future/runoff.ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ggarrange(plotlist = listgrob, ncol = 2, nrow = 2, labels = c("SSP126-2050","SSP370-2050","SSP126-2100","SSP370-2100"))
wr.sf <- readRDS("WR/wr.sf.rds")
ssp126.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp126.runoff.wr.mean <- ssp126.runoff.wr %>% mean()
ssp126.runoff.wr.mean.crop <- raster::crop(ssp126.runoff.wr.mean,c(0,35,55,73))
plot(ssp126.runoff.wr.mean.crop)
ssp126.runoff <- raster::extract(ssp126.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff) <- c("gid","runoff")
ssp126.runoff$Runoff <- ssp126.runoff$runoff *365*24*60*60
saveRDS(ssp126.runoff,"WR/ssp126.runoff.2050.rds")
ssp370.runoff.wr <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(313:552))
ssp370.runoff.wr.mean <- ssp370.runoff.wr %>% mean()
ssp370.runoff.wr.mean.crop <- raster::crop(ssp370.runoff.wr.mean,c(0,35,55,73))
plot(ssp370.runoff.wr.mean.crop)
ssp370.runoff <- raster::extract(ssp370.runoff.wr.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff) <- c("gid","runoff")
ssp370.runoff$Runoff <- ssp370.runoff$runoff *365*24*60*60
saveRDS(ssp370.runoff,"WR/ssp370.runoff.2050.rds")
wr.sf <- readRDS("WR/wr.sf.rds")
raster.ssp126.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp126_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp126.runoff.2100.mean <- raster.ssp126.runoff.2100 %>% mean()
raster.ssp126.runoff.2100.mean.crop <- raster::crop(raster.ssp126.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp126.runoff.2100.mean.crop)
ssp126.runoff.2100 <- raster::extract(raster.ssp126.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp126.runoff.2100) <- c("gid","runoff")
ssp126.runoff.2100$Runoff <- ssp126.runoff.2100$runoff *365*24*60*60
saveRDS(ssp126.runoff.2100,"WR/ssp126.runoff.2100.rds")
raster.ssp370.runoff.2100 <- raster::stack("CORDEX_future/mrros_Lmon_CNRM-CM6-1-HR_ssp370_r1i1p1f2_gr_201501-210012.nc", bands = c(805:1032))
raster.ssp370.runoff.2100.mean <- raster.ssp370.runoff.2100 %>% mean()
raster.ssp370.runoff.2100.mean.crop <- raster::crop(raster.ssp370.runoff.2100.mean,c(0,35,55,73))
plot(raster.ssp370.runoff.2100.mean.crop)
ssp370.runoff.2100 <- raster::extract(raster.ssp370.runoff.2100.mean,dplyr::select(wr.sf,c("gid","geom")),fun=mean,df = T, sp = T)
names(ssp370.runoff.2100) <- c("gid","runoff")
ssp370.runoff.2100$Runoff <- ssp370.runoff.2100$runoff *365*24*60*60
saveRDS(ssp370.runoff.2100,"WR/ssp370.runoff.2100.rds")
Future TNdep will largely depend on the public policies implemented under the scenarios SSP 1-2.6 and 3-7.0.
TNdep data are from 2000.
Emission reduction between 2005 and 2020 (European Environment Agency 2021):
Predictions for future emissions for NOx (summary for policymakers AR6 (Masson-Delmotte et al. 2021)):
Assumption for NH3:
ndep.wr <- readRDS("WR/ndep.wr.rds")
ndep.sum <- ndep.wr@data %>% apply(2, sum)
(ndep.sum["woxn"]+ndep.sum["doxn"])/ndep.sum["tndep"]
(ndep.sum["wrdn"]+ndep.sum["drdn"])/ndep.sum["tndep"]
# SSP126 2050
ssp126.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2050$oxn <- ndep.wr$woxn*0.55*0.8 + ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$rdn <- ndep.wr$wrdn*0.9*0.8 + ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2050
ssp126.ndep.2050$TNdep <- ssp126.ndep.2050$oxn + ssp126.ndep.2050$rdn
saveRDS(ssp126.ndep.2050,"WR/ssp126.ndep.2050.rds")
# SSP126 2100
ssp126.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp126.ndep.2100$oxn <- ndep.wr$woxn*0.55*0.8+ndep.wr$doxn*0.55*0.8 # -45% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$rdn <- ndep.wr$wrdn*0.9*0.8+ndep.wr$drdn*0.9*0.8 # -10% for decrease between 2000 and 2020, -20% between 2020 and 2100
ssp126.ndep.2100$TNdep <- ssp126.ndep.2100$oxn + ssp126.ndep.2100$rdn
saveRDS(ssp126.ndep.2100,"WR/ssp126.ndep.2100.rds")
# SSP370 2050
ssp370.ndep.2050 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2050$oxn <- ndep.wr$woxn*0.55*1.37+ndep.wr$doxn*0.55*1.37 # -45% for decrease between 2000 and 2020, 137% increase between 2020 and 2100
ssp370.ndep.2050$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2050$TNdep <- ssp370.ndep.2050$oxn + ssp370.ndep.2050$rdn
saveRDS(ssp370.ndep.2050,"WR/ssp370.ndep.2050.rds")
# SSP370 2100
ssp370.ndep.2100 <- data.frame(gid = ndep.wr$gid)
ssp370.ndep.2100$oxn <- ndep.wr$woxn*0.55*1.9+ndep.wr$doxn*0.55*1.9 # -45% for decrease between 2000 and 2020, 190% increase between 2020 and 2100
ssp370.ndep.2100$rdn <- ndep.wr$wrdn*0.9+ndep.wr$drdn*0.9 # -10% for decrease between 2000 and 2020, no change between 2020 and 2100
ssp370.ndep.2100$TNdep <- ssp370.ndep.2100$oxn + ssp370.ndep.2100$rdn
saveRDS(ssp370.ndep.2100,"WR/ssp370.ndep.2100.rds")
The datasets were gathered into a single dataframe.
ndvi.ssp126 <- readRDS("WR/ndvi.ssp126.rds")
ndvi.ssp370 <- readRDS("WR/ndvi.ssp370.rds")
ssp126.runoff.2050 <- readRDS("WR/ssp126.runoff.2050.rds")
ssp370.runoff.2050 <- readRDS("WR/ssp370.runoff.2050.rds")
ssp126.ndep.2050 <- readRDS("WR/ssp126.ndep.2050.rds")
ssp370.ndep.2050 <- readRDS("WR/ssp370.ndep.2050.rds")
ndvi.ssp126.2100 <- readRDS("WR/ndvi.ssp126.2100.rds")
ndvi.ssp370.2100 <- readRDS("WR/ndvi.ssp370.2100.rds")
ssp126.runoff.2100 <- readRDS("WR/ssp126.runoff.2100.rds")
ssp370.runoff.2100 <- readRDS("WR/ssp370.runoff.2100.rds")
ssp126.ndep.2100 <- readRDS("WR/ssp126.ndep.2100.rds")
ssp370.ndep.2100 <- readRDS("WR/ssp370.ndep.2100.rds")
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
# 2050 -----
wr.sf.126.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126, !c("Arable", "Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2050)), by = "gid") %>% merge(ssp126.ndep.2050)
wr.sf.370.2050 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2050)), by = "gid") %>% merge(ssp370.ndep.2050)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2050))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2050))==F)
wr.sf.126.2050 <- wr.sf.126.2050[-na.lines.126,]
wr.sf.370.2050 <- wr.sf.370.2050[-na.lines.370,]
wr.sf.126.2050$logRunoff <- wr.sf.126.2050$Runoff %>% log()
wr.sf.370.2050$logRunoff <- wr.sf.370.2050$Runoff %>% log()
saveRDS(wr.sf.126.2050, "WR/wr.sf.126.2050.rds")
saveRDS(wr.sf.370.2050, "WR/wr.sf.370.2050.rds")
# 2100 ----
wr.sf.126.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp126.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp126.runoff.2100)), by = "gid") %>% merge(ssp126.ndep.2100)
wr.sf.370.2100 <- dplyr::select(wr.sf.95, c("gid","Bog","Runoff","logRunoff","Temp","Precip","NDVI","Arable","TNdep")) %>% setNames(c("gid","Bog","Runoff95","logRunoff95","Temp95","Prec95","NDVI95","Arable","TNdep95","geometry")) %>% merge(select(ndvi.ssp370.2100, !c("Arable","Runoff")), by = "gid") %>% merge(st_drop_geometry(st_as_sf(ssp370.runoff.2100)), by = "gid") %>% merge(ssp370.ndep.2100)
na.lines.126 <- which(complete.cases(st_drop_geometry(wr.sf.126.2100))==F)
na.lines.370 <- which(complete.cases(st_drop_geometry(wr.sf.370.2100))==F)
wr.sf.126.2100 <- wr.sf.126.2100[-na.lines.126,]
wr.sf.370.2100 <- wr.sf.370.2100[-na.lines.370,]
wr.sf.126.2100$logRunoff <- wr.sf.126.2100$Runoff %>% log()
wr.sf.370.2100$logRunoff <- wr.sf.370.2100$Runoff %>% log()
saveRDS(wr.sf.126.2100, "WR/wr.sf.126.2100.rds")
saveRDS(wr.sf.370.2100, "WR/wr.sf.370.2100.rds")
Then we plot the difference between the forecasted variables in 2050 and 2100, compared to 1995.
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
t1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 1-2.6, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
t3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2050", x = "Temperature 1995 (°C)")+theme_minimal()
t4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Temp95, y = Temp))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Temperature SSP 3-7.0, 2100", x = "Temperature 1995 (°C)")+theme_minimal()
ggarrange(t1,t2,t3,t4, nrow = 2, ncol = 2)
p1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 1-2.6, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y= "Precipitations SSP 1-2.6, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2050", x = "Precipitations 1995 (mm/y)")+theme_minimal()
p4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = Prec95, y = Precip))+geom_abline(slope = 1, intercept = 0, col = "red")+
labs(y = "Precipitations SSP 3-7.0, 2100", x = "Precipitations 1995 (mm/y)")+theme_minimal()
ggarrange(p1,p2,p3,p4, nrow = 2, ncol = 2)
# d1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
# d4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = TNdep95, y = TNdep))+geom_abline(slope = 1, intercept = 0, col = "red")+theme_minimal()
#
# ggarrange(d1,d2,d3,d4, nrow = 2, ncol = 2)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
n1 <- ggplot(wr.sf.126.2050)+geom_point(aes(x = NDVI95, y = NDVI))+ geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2050")+ theme_minimal()
n2 <- ggplot(wr.sf.126.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 1-2.6, 2100")+theme_minimal()
n3 <- ggplot(wr.sf.370.2050)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2050")+theme_minimal()
n4 <- ggplot(wr.sf.370.2100)+geom_point(aes(x = NDVI95, y = NDVI))+geom_abline(slope = 1, intercept = 0, col = "red")+
ylim(0,1)+labs(y = "NDVI SSP 3-7.0, 2100")+theme_minimal()
ggarrange(n1,n2,n3,n4, nrow = 2, ncol = 2)
glm.ndvi <- readRDS("WR/glm.ndvi.rds")
fennoscandia <- readRDS("fennoscandia.33.rds")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
par(mfrow = c(2,2))
xtemp1 <- seq(min(wr.sf.126.2050$Temp), max(wr.sf.126.2050$Temp), 0.1)
yprecip1 <- seq(min(wr.sf.126.2050$logPrecip), max(wr.sf.126.2050$logPrecip), 0.1)
data.fit1 <- expand.grid(Temp = xtemp1, logPrecip = yprecip1)
mtrx3d1 <- predict(glm.ndvi, newdata = data.fit1, type = "response")
mtrx.melt1 <- cbind(data.fit1, mtrx3d1)
names(mtrx.melt1) <- c("Temp","logPrecip","NDVI")
fcp1 <- ggplot(mtrx.melt1,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2050")
xtemp2 <- seq(min(wr.sf.126.2100$Temp), max(wr.sf.126.2100$Temp), 0.1)
yprecip2 <- seq(min(wr.sf.126.2100$logPrecip), max(wr.sf.126.2100$logPrecip), 0.1)
data.fit2 <- expand.grid(Temp = xtemp2, logPrecip = yprecip2)
mtrx3d2 <- predict(glm.ndvi, newdata = data.fit2, type = "response")
mtrx.melt2 <- cbind(data.fit2, mtrx3d2)
names(mtrx.melt2) <- c("Temp","logPrecip","NDVI")
fcp2 <- ggplot(mtrx.melt2,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.126.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 1-2.6, 2100")
xtemp3 <- seq(min(wr.sf.370.2050$Temp), max(wr.sf.370.2050$Temp), 0.1)
yprecip3 <- seq(min(wr.sf.370.2050$logPrecip), max(wr.sf.370.2050$logPrecip), 0.1)
data.fit3 <- expand.grid(Temp = xtemp3, logPrecip = yprecip3)
mtrx3d3 <- predict(glm.ndvi, newdata = data.fit3, type = "response")
mtrx.melt3 <- cbind(data.fit3, mtrx3d3)
names(mtrx.melt3) <- c("Temp","logPrecip","NDVI")
fcp3 <- ggplot(mtrx.melt3,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2050, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2050")
xtemp4 <- seq(min(wr.sf.370.2100$Temp), max(wr.sf.370.2100$Temp), 0.1)
yprecip4 <- seq(min(wr.sf.370.2100$logPrecip), max(wr.sf.370.2100$logPrecip), 0.1)
data.fit4 <- expand.grid(Temp = xtemp4, logPrecip = yprecip4)
mtrx3d4 <- predict(glm.ndvi, newdata = data.fit4, type = "response")
mtrx.melt4 <- cbind(data.fit4, mtrx3d4)
names(mtrx.melt4) <- c("Temp","logPrecip","NDVI")
fcp4 <- ggplot(mtrx.melt4,aes(x=Temp, y = logPrecip, z = NDVI)) + stat_contour(geom="polygon", aes(fill = ..level..))+
geom_tile(aes(fill = NDVI))+
stat_contour(bins = 15)+
geom_point(data = wr.sf.370.2100, aes(x=Temp, y = logPrecip), size = 1)+
scale_fill_viridis_c()+
labs(title = "Contour plot NDVI SSP 3-7.0, 2100")
grid.arrange(fcp1,fcp2,fcp3,fcp4,nrow=2,ncol=2)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
mn0 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI95, col = NDVI95)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI 1995", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn0, filename = "WR/NDVI_1995.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn1 <- ggplot(wr.sf.126.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn1, filename = "WR/NDVI_SSP126_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn2 <- ggplot(wr.sf.126.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP126-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn2, filename = "WR/NDVI_SSP126_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn3 <- ggplot(wr.sf.370.2050)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2050", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn3, filename = "WR/NDVI_SSP370_2050.png",dpi = "retina", width = 12, height = 16, units = "cm")
mn4 <- ggplot(wr.sf.370.2100)+geom_sf(aes(fill = NDVI, col = NDVI)) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0.5, rev = T , aesthetics = c("fill","col"), name = "NDVI SSP370-2100", limits = c(0,1))+ theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = mn4, filename = "WR/NDVI_SSP370_2100.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/NDVI_1995.png","white.png","WR/NDVI_SSP126_2050.png", "WR/NDVI_SSP126_2100.png","WR/NDVI_SSP370_2050.png","WR/NDVI_SSP370_2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 3, align = "v", greedy = T)
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
diff.var <- function(sf, vars1, vars2){
sf2 <- sf
for (i in 1:length(vars2)){
newname <- paste("diff", vars2[i], sep = "")
x <- sf2[,vars2[i]] %>% st_drop_geometry()
y <- sf2[,vars1[i]] %>% st_drop_geometry()
z <- x-y
names(z) <- newname
sf2 <- cbind(sf2, z)
}
return(sf2)
}
wr.sf.126.2050 <- diff.var(wr.sf.126.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2050 <- diff.var(wr.sf.370.2050, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.126.2100 <- diff.var(wr.sf.126.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
wr.sf.370.2100 <- diff.var(wr.sf.370.2100, vars1 = c("Temp95","Prec95","NDVI95","Runoff95","TNdep95"), vars2 = c("Temp","Precip","NDVI","Runoff","TNdep"))
list.sf <- list(wr.sf.126.2050, wr.sf.126.2100, wr.sf.370.2050, wr.sf.370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
all.diff <- data.frame(gid = c(wr.sf.126.2050$gid, wr.sf.370.2050$gid, wr.sf.126.2100$gid, wr.sf.370.2100$gid))
for(x in c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")){
all.vec <- c()
for(i in 1:length(list.sf)){
vec <- select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
# x1 <- filter(wr.sf.126.2050, NDVI95 > 0.4)
# x2 <- filter(wr.sf.370.2050, NDVI95 > 0.4)
# x3 <- filter(wr.sf.126.2100, NDVI95 > 0.4)
# x4 <- filter(wr.sf.370.2100, NDVI95 > 0.4)
# all.ndvi <- c(x1$NDVI/x1$NDVI95*100-100,x2$NDVI/x2$NDVI95*100-100, x3$NDVI/x3$NDVI95*100-100, x4$NDVI/x4$NDVI95*100-100)
# palettes <- data.frame(diffTemp = brewer.pal(name = "RdYlBu", n = 9)[c(1,5,9)],
# diffPrecip = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffNDVI = brewer.pal(name = "RdYlGn", n = 9)[c(9,5,1)],
# diffRunoff = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)],
# diffTNdep = brewer.pal(name = "RdYlBu", n = 9)[c(9,5,1)])
legends <- c("Difference Temperature (°C)", "Difference Precipitation (mm/y)", "Difference NDVI", "Difference Runoff (mm/y)","DIfference Nitrogen Deposition (ug/m2)")
diffs <- c("diffTemp","diffPrecip","diffNDVI","diffRunoff","diffTNdep")
rev.pal <- c(T,F,T,F,T)
for(i in c(1:length(list.sf))){
for(x in diffs[3]){
g <- ggplot(list.sf[[i]])+geom_sf(aes(fill = .data[[x]], col = .data[[x]])) +
scale_colour_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = rev.pal[grep(x, diffs)] , aesthetics = c("fill","col"), name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
# scale_fill_gradient2(low=palettes[3, grep(x,names(palettes))],mid=palettes[2, grep(x,names(palettes))],high = palettes[1, grep(x,names(palettes))],
# midpoint = 0, aesthetics = c("fill","col"),
# name = paste(legends[grep(x, diffs)],"\n",names(list.sf)[i], sep = ""), limits = c(min.lim[[x]], max.lim[[x]])) +
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave(plot = g, filename = paste("WR/", names(list.sf)[i], ".", x, ".png",sep = ""),dpi = "retina", width = 12, height = 16, units = "cm")
}
}
white <- ggplot()+theme_void()
ggsave(plot = white, filename = "white.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/SSP126.2050.diffTemp.png", "WR/SSP126.2100.diffTemp.png","WR/SSP370.2050.diffTemp.png","WR/SSP370.2100.diffTemp.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffPrecip.png","WR/SSP126.2100.diffPrecip.png","WR/SSP370.2050.diffPrecip.png","WR/SSP370.2100.diffPrecip.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffNDVI.png","WR/SSP126.2100.diffNDVI.png","WR/SSP370.2050.diffNDVI.png","WR/SSP370.2100.diffNDVI.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffRunoff.png","WR/SSP126.2100.diffRunoff.png","WR/SSP370.2050.diffRunoff.png","WR/SSP370.2100.diffRunoff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
listgrob <- list("WR/SSP126.2050.diffTNdep.png","WR/SSP126.2100.diffTNdep.png","WR/SSP370.2050.diffTNdep.png","WR/SSP370.2100.diffTNdep.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
The average concentration of total organic carbon in the coastal drainage basins was predicted based on the SELM with 5 predictors, fitted on the Fennoscandian dataset. The difference between the predicted log(TOC) for 2050 and 2100 under the SSP 1-2.6 and 3-7.0 scenarios, compared to the log(TOC) predicted for 1995, was computed and plotted.
fennoscandia <- readRDS("fennoscandia.33.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2050.rds")
temperature <- c(fennoscandia$Temp, sem.pred.wr$Temp,wr.ssp126.2050$Temp,wr.ssp370.2050$Temp, wr.ssp126.2100$Temp,wr.ssp370.2100$Temp)
precipitation <- c(fennoscandia$Precip, sem.pred.wr$Precip,wr.ssp126.2050$Precip,wr.ssp370.2050$Precip, wr.ssp126.2100$Precip,wr.ssp370.2100$Precip)
ndvi <- c(fennoscandia$NDVI, sem.pred.wr$NDVI,wr.ssp126.2050$NDVI,wr.ssp370.2050$NDVI, wr.ssp126.2100$NDVI,wr.ssp370.2100$NDVI)
runoff <- c(fennoscandia$Runoff, sem.pred.wr$Runoff,wr.ssp126.2050$Runoff,wr.ssp370.2050$Runoff, wr.ssp126.2100$Runoff,wr.ssp370.2100$Runoff)
tndep <- c(fennoscandia$TNdep, sem.pred.wr$TNdep,wr.ssp126.2050$TNdep,wr.ssp370.2050$TNdep, wr.ssp126.2100$TNdep,wr.ssp370.2100$TNdep)
toc <- c(fennoscandia$TOC, sem.pred.wr$TOC95,wr.ssp126.2050$TOC,wr.ssp370.2050$TOC, wr.ssp126.2100$TOC,wr.ssp370.2100$TOC)
neworder <- c("Fennoscandia","WR1995","SSP126.2050","SSP126.2100","SSP370.2050","SSP370.2100")
model <- factor(c(rep("Fennoscandia",dim(fennoscandia)[1]),
rep("WR1995",dim(sem.pred.wr)[1]),
rep("SSP126.2050", dim(wr.ssp126.2050)[1]),
rep("SSP370.2050", dim(wr.ssp370.2050)[1]),
rep("SSP126.2100", dim(wr.ssp126.2100)[1]),
rep("SSP370.2100", dim(wr.ssp370.2100)[1])),
levels = neworder)
compare.ssp <- data.frame(value = c(temperature, precipitation, runoff, ndvi, toc), variable = c(rep("Temperature (°C)", length(temperature)), rep("Precipitation (mm)", length(precipitation)), rep("Runoff (mm/y)",length(runoff)), rep("NDVI",length(ndvi)), rep("TOC (mg/L)", length(toc))), model = rep(model,5))
ggplot(compare.ssp)+geom_boxplot(aes(y= value,col=variable))+
labs(y="", x = "")+scale_color_viridis_d(end = 0.8)+
theme_minimal()+
theme(legend.position = "", axis.text.x = element_blank())+
facet_grid(cols = vars(model), rows = vars(factor(variable, levels = c("Temperature (°C)","Precipitation (mm)","Runoff (mm/y)","NDVI","TOC (mg/L)"))), scales = "free_y", shrink = T)
ggsave("WR/compare.ssp.png", width = 30, height = 25, units = "cm", dpi = "retina")
knitr::include_graphics("WR/compare.ssp.png")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2050 <- readRDS("WR/wr.sf.126.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2050,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2050 <- readRDS("WR/wr.sf.370.2050.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2050,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2050,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2050.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.126.2100 <- readRDS("WR/wr.sf.126.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp126 <- predict(sem.fennoscandia.5,wr.sf.126.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp126 <- cbind(wr.sf.126.2100,wr.pred.ssp126)
wr.ssp126$TOC <- exp(wr.ssp126$fit)
wr.ssp126 <- merge(wr.ssp126, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp126,"WR/wr.ssp126.2100.rds")
sem.fennoscandia.5 <- readRDS("sem.fennoscandia.5.Rdata")
wr.sf.370.2100 <- readRDS("WR/wr.sf.370.2100.rds")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
fm <- logTOC~NDVI+logRunoff+Bog+Arable+TNdep
wr.pred.ssp370 <- predict(sem.fennoscandia.5,wr.sf.370.2100,wr.kmat, pred.type = "TS",legacy.mixed = F, zero.policy = T) %>% as.data.frame()
wr.ssp370 <- cbind(wr.sf.370.2100,wr.pred.ssp370)
wr.ssp370$TOC <- exp(wr.ssp370$fit)
wr.ssp370 <- merge(wr.ssp370, select(st_drop_geometry(sem.pred.wr),c("gid","wr.pred.fit","TOC95")))
saveRDS(wr.ssp370,"WR/wr.ssp370.2100.rds")
The average TOC concentration was plotted.
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds")
list.sf <- list(wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
for(i in names(list.sf)){
list.sf[[i]]$diffTOC <- list.sf[[i]]$TOC - list.sf[[i]]$TOC95
list.sf[[i]]$diff.log <- list.sf[[i]]$fit - list.sf[[i]]$wr.pred.fit
list.sf[[i]]$exp.diff.log <- exp(list.sf[[i]]$diff.log)
list.sf[[i]]$percent.change <- (list.sf[[i]]$exp.diff.log - 1) * 100
}
all.diff <- data.frame(gid = c(wr.ssp126.2050$gid, wr.ssp370.2050$gid, wr.ssp126.2100$gid, wr.ssp370.2100$gid))
for(x in c("fit","diffTOC","diff.log","exp.diff.log","percent.change")){
all.vec <- c()
for(i in names(list.sf)){
vec <- dplyr::select(st_drop_geometry(list.sf[[i]]),all_of(x))
all.vec <<- rbind(all.vec,vec)
}
all.diff <- cbind(all.diff, all.vec)
}
max.lim <- apply(all.diff, 2, max, na.rm = T)
min.lim <- apply(all.diff, 2, min, na.rm = T)
t0 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=wr.pred.fit, col= wr.pred.fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "logTOC in 1995")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/toc_1995.png", t0 ,dpi = "retina", width = 12, height = 16, units = "cm")
t1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "a) Forecast logTOC\n SSP1-2.6 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
# scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast logTOC\n SSP1-2.6 - 2050", mid = (max.lim[2]-min.lim[2])/2, limits = c(min.lim[2],max.lim[2]), rev = T) +
# theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2050.png", t1 ,dpi = "retina", width = 12, height = 16, units = "cm")
t2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "b) Forecast logTOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp126.2100.png", t2 ,dpi = "retina", width = 12, height = 16, units = "cm")
t3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "c) Forecast logTOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2050.png", t3 ,dpi = "retina", width = 12, height = 16, units = "cm")
t4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=fit, col= fit), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[2],max.lim[2]),
name = "d) Forecast logTOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ggsave("WR/logTOC_ssp370.2100.png", t4 ,dpi = "retina", width = 12, height = 16, units = "cm")
knitr::include_graphics("WR/toc_1995.png")
list.grob <- c("WR/logTOC_ssp126.2050.png", "WR/logTOC_ssp126.2100.png","WR/logTOC_ssp370.2050.png","WR/logTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
p1 <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "a) Forecasted percent increase in [TOC]\n SSP1-2.6 - 2050")+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] increase in %\n SSP1-2.6 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p2 <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] increase in %\n SSP1-2.6 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p3 <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] increase in %\n SSP3-7.0 - 2050", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
p4 <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] increase in %\n SSP3-7.0 - 2100", mid = 1, limits = c(min.lim[5],max.lim[5]), rev = T) +
# scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],max.lim[5]),
# name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
theme_minimal(base_size = 6)+
theme(legend.position = "bottom")
pp <- grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
ggsave("forecast.percent.change.toc.png", pp, width = 12, height = 16, units = "cm")
highchangeid <- subset(list.sf$SSP370.2100, exp.diff.log > 25)$gid
pp1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "a) Forecasted percent change TOC\n SSP1-2.6 - 2050")+
geom_sf(data = subset(list.sf[["SSP126.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "b) Forecasted percent change TOC\n SSP1-2.6 - 2100")+
geom_sf(data = subset(list.sf[["SSP126.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "c) Forecasted percent change TOC\n SSP3-7.0 - 2050")+
geom_sf(data = subset(list.sf[["SSP370.2050"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
pp4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log < 20), aes(fill=exp.diff.log, col= exp.diff.log), size= 0.05)+
scale_fill_distiller(type = "seq", palette = 8, aesthetics = c("colour","fill"), direction = 1, limits = c(min.lim[5],20),
name = "d) Forecasted percent change TOC\n SSP3-7.0 - 2100")+
geom_sf(data = subset(list.sf[["SSP370.2100"]], exp.diff.log > 20), aes(fill=exp.diff.log), col= "black", fill = "black", size= 0.05)+
theme_minimal(base_size = 10)+
theme(legend.position = "bottom")
ppp <- grid.arrange(pp1, pp2, pp3, pp4, ncol = 2, nrow = 2)
ggsave("WR/forecast.prop.change.png", ppp, width = 12, height = 16)
c1 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Forecast [TOC] change in %\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2050.png", c1 ,dpi = "retina", width = 12, height = 16, units = "cm")
c2 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Forecast [TOC] change in %\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP126.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp126.2100.png", c2 ,dpi = "retina", width = 12, height = 16, units = "cm")
c3 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Forecast [TOC] change in %\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2050"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2050.png", c3 ,dpi = "retina", width = 12, height = 16, units = "cm")
c4 <- ggplot()+
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change < 307), aes(fill=percent.change, col= percent.change), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Forecast [TOC] change in %\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[6],307), rev = T) +
geom_sf(data = subset(list.sf[["SSP370.2100"]], percent.change > 307), fill = "black", col = "black")+
theme_minimal(base_size = 12)+
theme(legend.position = "bottom")
ggsave("WR/percentTOC_ssp370.2100.png", c4 ,dpi = "retina", width = 12, height = 16, units = "cm")
list.grob <- c("WR/percentTOC_ssp126.2050.png", "WR/percentTOC_ssp126.2100.png","WR/percentTOC_ssp370.2050.png","WR/percentTOC_ssp370.2100.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
plot_grid(plotlist = list.grob, nrow = 2, ncol = 2, labels = c("a)","b)","c)","d)"), label_size = 10)
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="a) Difference logTOC\n SSP1-2.6 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="b) Difference logTOC\n SSP1-2.6 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="c) Difference logTOC\n SSP3-7.0 - 2050", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=diff.log, col= diff.log), size= 0.05)+
scale_fill_continuous_divergingx(palette = "RdYlBu", aesthetics = c("fill","col"),name="d) Difference logTOC\n SSP3-7.0 - 2100", mid = 0, limits = c(min.lim[4],max.lim[4]), rev = T) +
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.diff.log.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/SSP126.2050.diff.log.png", "WR/SSP126.2100.diff.log.png","WR/SSP370.2050.diff.log.png","WR/SSP370.2100.diff.log.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
ssp <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("ssp.toc.png",ssp, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(ssp)
for(i in names(list.sf)){
list.sf[[i]]$sign.diff <- sign(list.sf[[i]]$diffTOC) %>% as.factor()
}
g <- ggplot(list.sf[["SSP126.2050"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 1-2.6 2050")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2050.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP126.2100"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 1-2.6 2100")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP126.2100.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2050"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 3-7.0 2050")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2050.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
g <- ggplot(list.sf[["SSP370.2100"]])+
geom_sf(aes(fill=sign.diff, col= sign.diff), size= 0.05)+
scale_color_manual(values = c("lightblue","firebrick"), labels = c("Decreased [TOC]", "Increased [TOC]"), aesthetics = c("fill","col"), name = "SSP 3-7.0 2100")+
theme_minimal(base_size = 10)+theme(legend.position = "bottom")
ggsave(plot = g, filename = "WR/SSP370.2100.sign.diff.png",dpi = "retina", width = 12, height = 16, units = "cm")
listgrob <- list("WR/SSP126.2050.sign.diff.png", "WR/SSP126.2100.sign.diff.png","WR/SSP370.2050.sign.diff.png","WR/SSP370.2100.sign.diff.png") %>% lapply(readPNG) %>% lapply(rasterGrob)
sdf <- plot_grid(plotlist = listgrob, ncol= 2, nrow = 2, align = "v", greedy = T)
save_plot("sign.diff.toc.png",sdf, ncol = 2, nrow = 2, base_height = 4, base_width = 3)
print(sdf)
Most coastal drainage basins in Sweden and Finland drain into to the Baltic Sea (see Supplementary 2), while Norwegian sea watershed drain in Skagerrak, in the North sea, in the Norwegian Sea and in the Barent sea in the North (Sætre 2007). The drainage basin of the Baltic sea matches almost exactly the borders of Sweden and Finland : https://www.grida.no/resources/5324. All contribute to the Norwegian Coastal Current water masses. The estimated exported TOC was computed as:
\[ TOC_exp = [TOC] (mg.L^{-1}) \times Runoff (mm/y, or L.m^{-2}.y^{-1}) \times Watershed area (m^2)\]
Where the TOC concentration (past and future) was obtained by back-transforming the results of the model (log(TOC)) into TOC in mg/L. The result of the above equation, giving the amount of TOC exported in mg, is converted to Tg. Forecast are presented on Figure 19*.
wr.sf.95 <- readRDS("WR/wr.sf.95.rds")
wr.centroid <- st_centroid(wr.sf.95, of_largest_polygon = T) %>% dplyr::select("gid","geometry")
sem.pred.wr <- readRDS("WR/sem.pred.wr.rds")
nor <- st_read("Country_shapefile/norway.shp") %>% st_transform(st_crs(wr.centroid))
swe <- st_read("Country_shapefile/sweden.shp") %>% st_transform(st_crs(wr.centroid))
fin <- st_read("Country_shapefile/finland.shp") %>% st_transform(st_crs(wr.centroid))
wr.norge <- wr.centroid[nor,]
wr.norge$Nation <- "Norway"
wr.sweden <- wr.centroid[swe,]
wr.sweden$Nation <- "Sweden"
wr.finland <- wr.centroid[fin,]
wr.finland$Nation <- "Finland"
wr.nation <- rbind(wr.norge,wr.sweden,wr.finland)
wr.sf.95 <- wr.sf.95 %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid") %>% sp::merge(dplyr::select(st_drop_geometry(sem.pred.wr), c("gid","wr.pred.fit","TOC95")), by = "gid")
names(wr.sf.95[which(names(wr.sf.95) == "TOC95")]) <- "TOC"
wr.ssp126.2050 <- readRDS("WR/wr.ssp126.2050.rds")%>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2050 <- readRDS("WR/wr.ssp370.2050.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp126.2100 <- readRDS("WR/wr.ssp126.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
wr.ssp370.2100 <- readRDS("WR/wr.ssp370.2100.rds") %>% sp::merge(dplyr::select(st_drop_geometry(wr.nation), c("gid","Nation")), by = "gid")
list.sf <- list(wr.sf.95, wr.ssp126.2050, wr.ssp126.2100, wr.ssp370.2050, wr.ssp370.2100)
names(list.sf) <- c("Data1995","SSP126.2050", "SSP126.2100", "SSP370.2050", "SSP370.2100")
export <- function(sf){
discharge <- sf$Runoff*sf$area # Q in L/y
toc.export <- discharge * sf$TOC * 10^(-15) # amount of exported TOC in Tg/y
sf2 <- cbind(sf,discharge,toc.export)
return(sf2)
}
list.sf <- lapply(list.sf, export)
list.mean <- c()
for(i in 1:length(list.sf)){
list.mean[[i]] <- list.sf[[i]] %>% st_drop_geometry() %>% group_by(Nation) %>% summarise(total = sum(toc.export)) %>% as.data.frame()
}
names(list.mean) <- names(list.sf)
list.diff <- c()
for(i in 1:length(list.mean)){
df1 <- list.mean[[1]]
df2 <- list.mean[[i]]
diff <- df2[,2] - df1[,2]
x <- cbind(list.mean[[i]],diff)
names(x) <- c("Nation", names(list.mean)[i], paste("diff", names(list.mean)[i], sep = ""))
list.diff[[i]] <- x
}
toc.export.df <- list.diff[[1]] %>% merge(list.diff[[2]], by = "Nation") %>% merge(list.diff[[3]], by = "Nation") %>% merge(list.diff[[4]], by = "Nation") %>% merge(list.diff[[5]], by = "Nation") %>% melt(id.vars = c("Nation"), value.name = "TOC_Tg")
toc.export.df$Model <- NA
m.126 <- grep("126",toc.export.df$variable)
m.370 <- grep("370",toc.export.df$variable)
toc.export.df$Model[m.126] <- "SSP126"
toc.export.df$Model[m.370] <- "SSP370"
toc.export.df$Model[-c(m.126,m.370)] <- "Historical"
toc.export.df$Year <- NA
m.95 <- grep("95",toc.export.df$variable)
m.2050 <- grep("2050",toc.export.df$variable)
m.2100 <- grep("2100",toc.export.df$variable)
toc.export.df$Year[m.95] <- "1995"
toc.export.df$Year[m.2050] <- "2050"
toc.export.df$Year[m.2100] <- "2100"
toc.export.df$Diff <- NA
m.d <- grep("diff",toc.export.df$variable)
toc.export.df$Diff[m.d] <- "Difference"
toc.export.df$Diff[-m.d] <- "Absolute"
toc.export.sum <- toc.export.df %>% group_by(Model,Year,Diff) %>% summarise(TOC_Tg = sum(TOC_Tg))
toc.export.sum$Nation <- "Total"
toc.export.df <- rbind(toc.export.sum,toc.export.df)
toc.export.df <- transform(toc.export.df, Nation = factor(Nation, levels = c("Norway","Sweden","Finland","Total")))
toc.export.df$variable <- NULL
toc.export.df$TOC_Tg <- as.numeric(toc.export.df$TOC_Tg)
knitr::kable(toc.export.df) %>% kable_styling()
saveRDS(toc.export.df,"WR/toc.export.df.rds")
ggplot() +
geom_point(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, shape = Model)) +
geom_path(data = toc.export.df, aes(x=Year, y = TOC_Tg, col = Model, group = Model), size = 0.5) +
geom_abline(slope = 0, intercept = 0, col = "firebrick4", size = 0.3)+
scale_color_manual(values = c("firebrick4","steelblue4","darkorange"))+
geom_label(data = filter(toc.export.df, Model == "SSP126" & Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP126"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = -0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "SSP370"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Absolute"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.5, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
geom_label(data = filter(toc.export.df, Model == "Historical"& Diff == "Difference"),
aes(x=Year, y = TOC_Tg, col = Model, label = round(TOC_Tg,1)), nudge_y = 0.1, size = 2, label.padding = unit(0.2,"lines"), show.legend = F) +
labs(x="Year", y= "TOC (Tg)")+
facet_grid(cols=vars(Nation), rows = vars(Diff), scales = "free_y")+
theme_bw(base_size = 12)+
theme(legend.position = "bottom")
ggsave(filename = "WR/TOC_export.png", dpi = "retina", width = 20, height = 15, units = "cm")
toc.export.df <- readRDS("WR/toc.export.df.rds") %>% dplyr::select(c("Nation", "Model", "Year","Diff","TOC_Tg"))
toc.export.df %>% arrange(Model) %>%
knitr::kable(digits = 3,
caption = "TOC export by country",
col.names = c("Country","Model","Year","Absolute/differential export","TOC export (Tg)")) %>%
kable_styling(bootstrap_options = "bordered",
position = "center",
full_width = F)
knitr::include_graphics("WR/TOC_export.png")